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Long-term evolution of accretion discs in Be/X-ray binaries 
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ABSTRACT 

We numerically study the long-term evolution of the accretion disc around the 
neutron star in a coplanar Be/X-ray binary with a short period and a moderate 
eccentricity. From three dimensional Smoothed Particle Hydrodynamics simulations, 
we find that the disc evolves through three distinct phases, each characterized by 
different mass accretion patterns. In the first "developing phase", the disc is formed 
and develops towards a nearly Keplerian disc. It has a relatively large, double-peaked 
mass-accretion rate with the higher peak by the direct accretion at periastron, which is 
followed by the lower peak by the accretion induced by a one-armed spiral wave. In the 
second "transition phase", the disc is approximately Keplerian and grows with time. 
The mass-accretion rate increases as the disc grows. In the second phase, there is a 
transition in the mass accretion rate from a double peaked to a single peaked pattern. 
In the final quasi-steady state, the mass-accretion rate is on average balanced with 
the mass-transfer rate from the Be disc and exhibits a regular orbital modulation. In 
the quasi-steady state, the mass-accretion rate has a single peak by the wave-induced 
accretion as in a later stage of the transition phase. The orbital modulation of X-ray 
maxima could provide not only a circumstantial evidence for the persistent disc but 
also an observational diagnosis of the disc evolutionary state. 

Key words: accretion, accretion discs - hydrodynamics - methods: numerical - 
binaries: general - stars: emission-line, Be - X-rays: binaries 



1 INTRODUCTION 

The Be/X-ray binaries comprise the largest subclass of high- 
mass X-ray binaries (HMXBs). About two-thirds of the iden- 
tified systems fall into this category. These systems con- 
sist of, generally, a neutron star and a Be star with a cool 
(~ 10 4 A") equatorial disc, which is geometrically thin and 
nearly Keplerian. These systems are distributed over a wide 
range of orbital periods (10 d < P or b < 300 d) and eccentric- 
ities (e < 0.9). 

Most of the Be/X-ray binaries show only transient ac- 
tivity in the X-ray emission and are termed Be/X-ray tran- 
sients. Be/X-ray transients show periodical (Type I) out- 
bursts, which are separated by the orbital period and have 
the luminosity of Lx = 10 36_37 ergs _1 , and giant (Type 
II) outbursts of Lx > 10 38 ergs _1 with no orbital modula- 
tion. These X-ray outbursts are basically considered to re- 
sult from the accretion onto the neutron star from the Be 
disc. 

Recently, lHavasaki fc Okazakil j2004l) (hereafter, Pa- 
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per I) studied the accretion flow around the neutron 
star in Be/X-ray binaries, performing three-dimensional 
(3D) smoothed particle hydrodynamics (S PH) simulations 
llBenzl ll99Ct iBate. Bonnell fc Pried 1 19951) . in the frame- 
work of t he interaction betwe en the Be disc and the neu- 
tron star (Oka zaki et aljl2002l) . They found that the time- 
dependent accretion d isc is formed around the neutron star. 
lHavasaki fc Okazaki (I2005T) (hereafter, Paper II) further 
showed that the disc has a one-armed spiral structure in- 
duced by a phase-dependent mass transfer from the Be disc. 
The one-armed mode is excited and then damped every or- 
bit, which make an enhancement in the X-ray luminosity. 
These are, however, the results from simulations run over a 
period shorter than the viscous time-scale of the disc, and 
it was not known how the accretion disc would evolve over 
a period longer than the viscous time-scale. 



In this paper, we study the long-term evolution of the 
accretion disc around the neutron star in Be/X-ray binaries, 
using a 3D SPH code. The simulations are performed using 
the basically the same method as in paper I. We describe 
our numerical results in Sections 2 and 3, and discuss the 
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Table 1. Summary of model simulations. The initial number of particles is 301 and the mean mass injection rate is ~ 2.5 X 
lO -11 p-iiM0yr — 1 in all simulations, where p_ n is the highest local density in the Be-star disc normalized by 10 _11 g, cm 3 . The 
first column represents the model. The second column is the polytropic exponent Y and the third column is the number of SPH particles 
at the end of the run. In the fourth column, the inner radius of the simulation region is described. The last column is the mean accretion 
rate for the last one orbital period. 



Model 


Polytropic exponent 


A^PH 


Inner boundary 


Afacc 




r 


(final) 


r- m / a 


(p-iiMoyr" 1 ) 


1 


1.2 


66676 


3.0 x 10" 3 


1.7 X lO" 11 


2 


5/3 (adiabatic) 


22122 


3.0 x 10" 3 


2.5 X 10~ n 


3 


1 (isothermal) 


49950 


1.0 x 10" 2 


1.5 x lO^ 11 



disc evolution and its observational implications in Section 
4. Section 5 summarizes our conclusions. 



2 DISC STRUCTURE 

We have carried out simulations by using the same com- 
puter code as in papers I and II. In our code, the neutron 
star is modeled by a sink particle with a fixed accretion 
radius n n , while the accretion disc is modeled by an en- 
semble of non-selfgravitating gas particles. As in papers I 
and II, our simulations are for a coplanar system with a 
short period (P or b = 24.3 d) and a moderate eccentricity 
(e = 0.34), the parameters of which are for 4U 0115+63, one 
of the best studied Be/X-ray binaries. The phase-dependent, 
mass-transfer rate from the Be disc is emulated by injecting 
gas particles, based on the result from a high resolution SPH 
simul ation for the same binary configuration (Okaz aki et alJ 
2002). The SPH particles are assumed to have the initial 
temperature of half the effective temperature of the Be star 
and follow a polytropic equation of state specified by the 
exponent P. 

We consider the case P = 1.2 (model 1) as well as the 
adiabatic case P = 5/3 (model 2) and the isothermal case 
r = 1 (model 3) as two extreme cases. Models 1-3 in this pa- 
per correspond to models 1, 3 and 4 in paper I, respectively. 
The injection rate of the SPH particles in each model is 
three times as low as in corresponding simulation in paper I 
in order to make the computing time bearable. In models 1 
and 2, we adopt r ln = 3.0 x 10 _3 a, where a is the semi-major 
axis of the binary. In model 3, however, a larger accretion 
radius of r- m = 1-0 x 10 _2 a is adopted, in order to reduce an 
artificial effect of the inner simulation boundary. All models 
have the Shakura-Sunyaev viscosity parameter ass =0.1 
and are run for 45P rt>, which is twice, 10 times and 5 times 
as long as the viscous time-scale at the disc outer radius in 
models 1, 2 and 3, respectively (see Section 3.3). The num- 
ber of neighbor particles is always kept to about 50 particles 
throughout our simulations Details of the models are given 
in Table 1. Throughout this paper, the unit of time is P or b, 
unless noted otherwise. 

We perform three dimensional (3D) simulations, despite 
that they have less radial resolution than two dimensional 
(2D) simulations do with the same particle number. 3D cal- 
culations have an advantage even for a geometrically thin 
disc in a coplanar binary. They enable us to evaluate more 
accurately the effect of the one-armed wave excited by the 
impulsive mass transfer from the Be disc than 2D calcula- 
tions do. As we will see in later sections, the one-armed wave 
plays an important role in the long-term evolution of the ac- 
cretion disc. 2D simulations would have artificially enhanced 



wave amplitud e because of the lack of vertical structure (e.g., 
lQgilvieH200ll) 1. 

In this section, we show the azimuthally averaged struc- 
tures, the non-axisymmetric structures and the vertical 
structures of the fully developed discs during the last one 
orbital period 44 < t < 45. 



2.1 Azimuthally averaged disc structure 

Fig.Qshows the radial structure of the accretion disc around 
the neutron star at apastron (t — 44.5). In each panel, the 
thick solid line, the dashed line and the thin solid line denote 
the surface density E, the azimuthal velocity normalized 
by the Keplarian velocity at the inner boundary and the ra- 
dial Mach number v r /c s , respectively. In the figure, the ve- 
locity components are averaged vertically and azimuthally, 
while the density is integrated vertically and averaged az- 
imuthally, where p_n is highest local density in the Be-star 
disc normalized by 10 _11 gcm -3 , a typical value for Be stars. 
It is noted from Fig.^that the disc is nearly Keplerian and 
the radial velocity component of the disc is highly subsonic 
in each model. It is also noted that the surface density in 
model 2 (adiabatic case) is the lowest because of the highest 
viscous stress due to the highest pressure. 

Fig. ^denotes the radial structures of the disc tempera- 
ture , the ratio of the smoothing length to the disc thickness 
and the ratio of the disc thickness to the radius at t = 44.5 in 
models 1-3. In each panel, the solid line, the dashed-dotted 
line and the dashed line denote the disc temperature T nor- 
malized by To = 1.3 x 10 4 K, which is the initial temperature 
of injected particles, the ratio of the smoothing length to the 
disc thickness h/H and the ratio of the disc thickness to the 
radius H/r, respectively. Here the disc thickness H is de- 
termined from the vertical distribution of the SPH particles 
by comparison with the analytical density distribution of 
a geometrically thin disc. The vertical disc structure and 
the procedure to determine H from the particle distribution 
will be given in Section [2.31 It is noted from the figure that 
the disc is fully resolved in the radial direction and is also 
resolved in the vertical direction except in the innermost 
region of the disc. 

In our simulations, the presence of the inner simulation 
boundary causes an artificial decrease in the density near the 
boundary, because the gas particles which pass through the 
boundary are removed from the simulation. As a result, the 
density distribution has a break near the inner boundary, 
as seen in Fig. We have found that our simulations are 
reliable for > 2n n . 
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Figure 1. Azimuthally averaged radial structures of the accretion disc at t=44.5 in models 1,2 and 3. In each panel, the thick solid 
line, the dashed line and the thin solid line denote the surface density S in units of p_ngcm 2 , the azimuthal velocity normalized by 
the Keplerian velocity at the inner boundary and the radial Mach number v r /c s , respectively. For comparison, the Keplerian velocity 
distribution is shown by the dash-dotted line. 
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Figure 2. Radial distributions of some characteristic quantities at t=44.5 in models 1,2 and 3. In each panel, the solid line, the dashed- 
dotted line and the dash line denote the disc temperature T normalized by To = 1.3 X 10 4 K which is the initial temperature of injected 
particles, the ratio of the smoothing length to the disc thickness h/H and the ratio of the disc thickness to the radius H/r, respectively. 



2.2 Non-axisymmetric disc structure 

It was found by paper II that an accretion disc around 
the neutron star in Be/X-ray binaries has a one-armed spi- 
ral structure induced by a phase-dependent mass transfer 
from the Be disc to the neutron star. In this subsection, we 
show the non-axisymmetric structure of the accretion disc 
for models 1-3. 

Fig IH gi ves snapshots of the accretion disc around the 
neutron star at t — 44.23 in model 1, at t = 44.12 in model 2 
and at t — 44.43 in model 3, respectively. Each time t is cho- 
sen such that the strength of the one-armed mode takes a 
maximum over the period of 44 < t < 45. The left pan- 
els show the contour maps of the surface density, whereas 
the non-axisymmetric components of the surface density and 
the velocity field are shown in the right panels. Annotated 
in each left panel are the time in units of P or b and the am- 
plitude of the one-armed wave Si, details of which will be 
described later. It is noted from the figure that the disc has 
a one-armed spiral structure regardless of the simulation pa- 
rameters, even after the disc fully evolves. The figure also 
suggests that the strength of one-armed mode correlates in- 
versely with the surface density. 



2.3 Vertical disc structure 

In this subsection, we compare the vertical distribution of 
SPH particles with the analytical density distribution of a 
geometrically thin disc. The vertical structure of such a disk 
is characterized by the disc scale-height (Gamma — 1) or 
thickness (Gammanel). In an isothermal disc, the vertical 
scale-height is analytically given by 

ffa = C S /Q. K , (1) 

whereas the disc thickness with the polytropic equation of 
state is defined by 

#a = ^ 7 " 1)/2 ((7 ~ l)/2 7 A' p )- 1/2 ^\ (2) 

where c s is the sound speed, Ox is the Keplerian frequency, 
po is the central density of the disc and K p is the polytropic 
constant. 

Fig. [3]shows the SPH particle distributions in the verti- 
cal direction in models 1-3. In each panel, the solid line and 
the dashed line denote the number of SPH particles normal- 
ized by the number of particles JV max at the disc mid-plane 
and the corresponding analytical distributions with given 
by equations Q or 

It is noted from the figure that the distributions of SPH 
particles are slightly broader than those of analytical distri- 
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Figure 3. Particle distributions in the vertical direction at t=44.5 in models 1-3. In each panel, the solid line and the dotted line 
denote the number of SPH particles normalized by the number of particles Nmax at the disc mid-plane and the corresponding analytical 
distribution, respectively. 



butions except for model 3 (isothermal disc). The ratio of the 
disc thickness obtained by the SPH particle distribution to 
the disc scale-hight or thickness by analytical ones H/H^. is 
1.12 in model 1, 1.32 in model 2 and 1.03 in model 3. The rea- 
son why the isothermal model has a better fit than models 
with Gammanel do is that the scale-height of the isother- 
mal disc depends only on the temperature, while models 
with Gammanel depend on the density, which has a numer- 
ical error. In any case, discs in our simulations are vertically 
resolved except in the innermost region. 



3 DISC EVOLUTION 

By analysing the multi- wave length long term monitor ing ob- 
servations of 4 U0115+63, iNeeueruela et alJ 11200 II) found 
that the Be star shows a quasi-cyclic activity (~ 3 — 5yr) 
due to the loss and reformation of its circumsteller disc. 
When the Be disc is lost, the mass-supply to the neutron 
star also has to halt. Since we have little knowledge on the 
mechanism that triggers the Be-disc loss episode, we have 
restricted our simulations for t < 3 yr ~ 45 P or b . 

In the first two subsections below, we will describe the 
evolution of the mass-accretion rate and the strength of 
modes for models 1-3. Then, we will show the evolution of 
the mass, radius and viscous time-scale of the accretion disc 
for these models. 



3.1 Mass accretion rate 

Fig-Elshows the evolution of the mass-accretion rate and the 
corresponding X-ray luminosity in models 1-3 for < t < 45. 
In each panel, the solid line and the dashed line denote the 
mass-accretion rate M aC c and the mean mass-transfer rate 
Mr , respectively, where Mr is ~ 2.5 x lO~ 11 p-nM j/r _1 . 
M acc in models 1, 2 and 3 are shown in the top panel, the 
middle panel and the bottom panel, respectively. 

In model 1 (top panel), the mass accretion rate has 
double peaks per orbit at an initial developing phase < 
t < 10. While the first peak is due to the direct accretion 
of particles with the low specific angular momentum, the 
second one is mainly caused by an inward propagation of 
m — 1 mode. After the disc is fully developed (i> 10), the 
mass accretion rate has a single peak per orbit only due 



to the wave induced accretion. In the long-term, the mass 
accretion rate gradually increases with time after t — 10 

In model 2 (middle panel), the mass accretion rate ba- 
sically shows the same behaviour as that of model 1 in the 
short term, except that its peak value is higher than the 
mean mass-transfer rate after the disc is developed. In the 
long term, the mass accretion rate gradually increases with 
time and saturates for t>20 with an orbital modulation. 
This is because the mass of gas transferred from the Be disc 
over one orbital period is balanced with that of the accret- 
ing gas on to the neutron star for t > 20. For example, as 
shown in Table Q the mean mass-accretion rate for the last 
one orbital period is 2.5 x 10~ p-uM^yr - , which approx- 
imately equals to the mean mass-transfer rate from the Be 
disc. Thus, we noted that the disc gets to the quasi-steady 
state for t > 20 in model 2. 

In model 3 (bottom panel), the mass accretion rate has 
double peaks per orbit throughout the simulation. Although 
the first peak of the mass accretion rate is remarkable for 
a developing phase of disc formation, it decrease with time 
after the disc is developed (t>20). The value of the first 
peak could include an artificial effect since it is related to 
the size of the inner simulation boundary ri„ — 1.0 x 10 _2 a. 

3.2 Mode strength 

As described in Section 12.21 the one-armed spiral wave is 
excited by the ram pressure due to phase-dependent mass- 
transfer from the Be disc regardless of the simulation param- 
eters. The strength of the mode can be defined in various 
ways by decomposing the surface density distribution into 
Fourier components which vary as exp (imtj)) with m being 
the azimutha l harm onic number. In paper II, we defined it 
after iLubowl (ll99lT) by integrating each Fourier component 
of the surface density over the whole disc and then calcu- 
lating its amplitude [see eq. (2) of paper II]. If the spiral 
wave is, however, tightly winding, this method will give sig- 
nificantly underestimated amplitudes. To avoid it, we define 
here the strength S m of the mode m as follows. 

S m (t) = -^-^dr(S c 2 os , m (r,t) + S 2 n , m (r, t)) 1/2 , (3) 

where S cos ,m and S s i n ,m are the azimuthal Fourier compo- 
nents of the surface density given by 
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Figure 5. Evolution of the mass-accretion rate M acc for models 1-3 with different polytropic exponents. In each panel, the solid line 
and the dashed horizontal line denote the mass-accretion rate and the averaged mass-transfer rate from the Be disc, respectively. The 
mass-accretion rate is measured in units of p_hMq yr — 1 . The right axis shows the X-ray luminosity corresponding to the mass-accretion 
rate with the X-ray emission efficiency T] = 1, where r) is defined by Lx = i)MxM a c C /fix. 



Sf, m (r,t) 



(1 + <5 m , ) J 



rdcf>T,(r,(f>,t)f(m<t>), 



(4) 



where / is either sin or cos function, and Md is the total 
disc mass given by 



M d = dr 



rd(f>T,(r, <f>, t). 



Fig. [(J shows the evolution of the strengths of several 
modes in models 1-3. In each panel, the solid thick line, the 
dotted line and the solid thin line show the strengths of 
m = 1, m = 2 and m = 3 modes, respectively. It is noted 
from the figure that the one-armed mode is always dominant 



among modes excited in the disc. In the short term, the 
m = 1 mode is amplified after every periastron passage by 
the ram pressure due to the phase-dependent mass transfer 
from the Be disc. The mode significantly decays by the next 
periastron passage. In the long-term, the amplitude of m = 
1 mode keeps decreasing until the disc reaches the quasi- 
steady state. 



3.3 Disc mass, radius and viscous time-scale 

Fig. I7f al shows the evolution of the disc mass in models 1-3 
over the period of < t < 45. In the figure, the thick solid 
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Figure 6. Evolution of several non-axisymmetric modes for models 1-3 with different polytropic exponent over the period of < t < 45. 
In each panel, the thick solid line, the dashed line and thin solid line denote the strengths of the m = 1, m = 2 and m = 3 modes, 
respectively. 



line, the dotted line and the thin solid line denote the disc 
mass in models 1, 2 and 3,respectively. In the short term, 
it modulates with the orbital phase. The disc mass starts 
increasing, when the mass-supply from the Be disc begins 
slightly before every periastron passage. It shows a rapid 
increase for a short while, followed by a gradual decrease 
by accretion with no mass-supply from the Be disc, which 
lasts until the next periastron passage. In the long term, 
the disc mass increases during the transition phase, because 
the averaged mass accretion rate on to the neutron star is 



smaller than the averaged mass transfer rate from the Be 
star. After the disc reaches a quasi-steady state, in which 
the mass accretion rate is, on average, balanced with the 
mass transfer rate, shows only a regular orbital modulation. 
For example, in model 2, the disc mass saturates at ~ 0.5 x 
1O" 12 M . 

In Fig. 0b) , we show the evolution of the disc radius 
over < t < 45 for models 1-3. As in paper I, the disc radius 
is measured by applying the following simple fitting function 
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Figure 4. Snapshots of the accretion disc for models 1-3. In each 
figure, the left panels show the surface density in a range of five 
orders of magnitude in the logarithmic scale, while the right pan- 
els show the non-axisymmetric components of the surface density 
(gray-scale plot) and the velocity field (arrows) in the linear scale. 
In the right panels, the region in gray (white) denotes the region 
with positive (negative) density enhancement. The periastron is 
in the ir-direction and the disc rotates counterclockwise. Anno- 
tated in each left panel are the time in units of P or b and the 
mode strength S\. 



E oc 



1 + (r/r d )" 



(5) 



to the radial distribution of the azimuthally-averaged sur- 
face density E. Here p and q are constants and r<j is the 
radius of accretion disc. In the early stage of disc formation, 
the disc is highly eccentric, for which the radius calcu- 
lated by equation JHJ is larger than that for a circular disc 
with the same disc mass. Since the disc eccentricity decreases 
with time during the developing phase, because of the vis- 
cous relaxation, so does the disc radius. After the disc is 
circularized and enters the transition phase, the disc radius 
increases with time until the disc reaches the quasi-steady 
state. 

Fig. |7Tc) shows the evolution of the viscous time-scale 
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Figure 8. Orbital-phase dependence of the peak of mass- 
accretion rate (the lower panel) and its frequency distribution 
(the upper panel). In the lower panel, the circles, the crosses and 
the triangles denote the orbital-phase dependence of the X-ray 
maxima for models 1-3, respectively. In the upper panel, the his- 
tograms filled by the hatched area, the solid area and the cross- 
hatched area show the frequency distributions corresponding to 
the X-ray maxima for models 1-3, respectively. 



of the disc over < t < 45 for models 1-3. The viscous 
time-scale r v i sc is estimated at the radius r/a = 0.06 in 
all models using azimuthally averaged quantities. Note that 
the smaller the polytropic exponent, the longer the viscous 
time-scale. Thus, the viscous time-scale in model 3 (r = 1) 
is the longest among three models studied in this paper. 
For example, at t — 45, r V i SC ~ 11.0, 5.5 and 24.0P O rb in 
models 1, 2 and 3, respectively. The disc evolution changes 
its trend dramatically before and after the viscous time-scale 
is reached, as shown in Figs. and and in Panels (b) and 
(c) of Fig.d 



4 DISCUSSION 

We have performed the three dimensional simulations of the 
long-term evolution of accretion disc around the neutron 
star in Be/X-ray binaries, taking into account the phase- 
dependent mass transfer from the Be disc. We have adopted 

the mass-transfer rate from a high-resolution simulation by 

II i I ° J 

Okazaki ct al. (2002) for a coplanar system with a short pe- 
riod (Porb = 24.3d) and moderate eccentricity (e = 0.34), 
which targeted 4U 01154-63, one of the best studied Be/X- 
ray binaries. 
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Figure 7. Evolution of (a) the disc mass in units of p_nl0 — 12 Mq, (b) the disc size normalized by the semi-major axis a and 
(c) the viscous time scale of the disc t„; 3c in units of the orbital period P or b for models 1-3 with different polytropic exponents over 
the period of < t < 45. In panel (c), each viscous time-scale is estimated by using t„; sc ~ r \l v ^ where v = ac s H is defined using 
azimuthally averaged quantities. 



Table 2. Characteristics of three evolutionary states of the accretion disc in Be/X-ray binaries. 'Wave-dependent phase' in the third 
column denotes that the peak phase of the wave-induced accretion depends on the excitation and propagation time-scales of the one-armed 
wave. 



Evolutionary state 


Accretion mechanism 


Orbital phase of a single or 
double accretion rate peaks 


Developing phase 
Transition phase 
Quasi-steady state 


Direct/Wavc-induced 
Direct (only in an earlier stage)/Wave-induced 
Wave-induced 


Periastron/Wavc-dcpcndcnt phase 
Pcriastron/ Wave-dependent phase 
Wave-dependent phase 
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We find that the disc evolves via three distinct phases, 
which consists of the first developing phase in which the 
disc is initially formed and developing, the second transition 
phase in which the disc is fully developed but still grows with 
time, and the final quasi-steady state in which the mass- 
transfer rate from the Be disc is balanced with the mass- 
accretion rate onto the neutron star. 

Throughout the developing phase and in an earlier stage 
of the transition phase, the mass-accretion rate has double 
peaks per orbit. The first peak, which is due to the direct 
accretion of gas particles with low specific angular momen- 
tum, occurs soon after periastron passage, whereas the sec- 
ond peak is due to the wave induced accretion and occurs 
at a phase dependent on the time-scales of excitation and 
propagation of the m — 1 wave. 

The first peak is dominant in the developing phase. But 
in the earlier stage of the transition phase, the predominance 
of the first peak decreases with time, and finally the second 
peak becomes stronger. In a later stage of the transition 
phase and in the quasi-steady state, the mass-accretion rate 
has a single peak due to the wave induced accretion. These 
characteristics are summarized in Table [5] 

Fig.|H]shows the peak phases of the mass accretion rate 
and their frequency distributions in models 1-3. Note that 
the accretion rate peaks are distributed over a wide range of 
orbital phase, depending on the polytropic exponent as well 
as the evolutionary state: In model 3 (T = 1), in which the 
direct accretion of particles of low specific angular momen- 
tum is the major accretion mechanism, the accretion rate 
peaks are concentrated around the periastron. In contrast, 
in model 1 (r = 1.2), in which the accretion is mainly by the 
m = 1 wave, the accretion rate peaks are distributed around 
the apastron with a large scatter. In model 2 (F = 5/3), 
in which the accretion is also due to the wave-induced ac- 
cretion, the peak phases are concentrated around 0.3P or b- 
The difference in the phases of accretion rate peaks between 
models 1, 2 and 3 results from the difference in the ma- 
jor accretion mechanism, whereas the peak-phase difference 
between models 1 and 2 is due to the difference between 
time-scales of excitation and propagation of the one-armed 
wave. 

The time-scale of wave propagation is roughly esti- 
mated by using the dispersion relati on of the o ne-armed 
oscillation in nearly Keplerian discs jKatolll983l) . Its fre- 
quency uj is written by u ~ — Q,K(k r H) 2 /2, where k r is 
the radial wave number. Then, the time-scale for a one- 
armed perturbation to travel across the disc is written by 
t w ~ r/\d,Lu/dk\ ~ assTVisc/27r. Thus, the time-scale of wave 
propagation r w is obtained as ~ 0.18P or b at t = 44.23 for 
model 1, ~ 0.08P O rb at t = 44.12 for model 2 and ~ 0.38P or b 
at t — 44.42 for model 3. Therefore, combing the time-scales 
for wave excitation and propagation, we expect that the X- 
ray maxima lag behind that of the mass-transfer rate by 
~ 0.41P or b, ~ 0.2P or b and ~ 0.8P or b for models 1-3, respec- 
tively. These results are in good agreement with the phase of 
the peak due to wave induced accretion (see Figs. and 01. 

The evolutionary time-scale of the disc depends on 
its viscous time-scale. The viscous time-scale of the 
disc increases with disc radius r^, because T v i S c(7"d) ~ 
assM-ffMA-dr^KM -1 oc rH 2 ' with (3 < being the 
temperature gradient in the radial direction. Thus, the larger 
the accretion disc, the longer the evolutionary time-scale. 



Given that the size of the accretion disc basically depends 
on the specific angular momentum of material transferred 
from the Be disc at periastron, which is larger for longer 
orbital period and/or smaller eccentricity, we expect that 
the accretion disc in Be/X-ray binaries evolves most slowly 
in systems with long period and low eccentricity, whereas 
it evolves most rapidly in short-period, highly eccentric sys- 
tems. 

On the other hand, the long-term mass supply to the 
neutron star in Be/X-ray binar ies is controlled by the q uasi- 
cycle of Be-star activity (see iNeeueruela et alJ 120011) for 
3 — 5 yr quasi-cycles of Be-disc formation and decay in 
4U 0115+63). Therefore, it is likely that accretion discs in 
systems with long period and low eccentricity show only the 
early evolutionary state, i.e., the developing phase and an 
early stage of the transition phase, whereas those in highly 
eccentric, short-period systems exhibit a variety of evolu- 
tionary states including the quasi-steady state. The observed 
X-ray behaviour can be used to probe the evolutionary state 
of the accretion disc in Be/X-ray binaries. 

Circinus X-l, which is a possible HMXB, has shown an 
interesting distribution of the X-ray maxima, a lthoug h the 
companion is unidentified yet. IClarkson et alJ J2004) have 
examined the distribution of the X-ray maxima in Circinus 
X-l by dividing mainly into a couple of different states ac- 
cording to the day numbers. In the high state, the X-ray 
maxima is distributed over a wide range of orbital phase 
with a peak at the orbital phase ~ 0.1P or b, whereas the 
distribution of the X-ray maxima is much more stable in 
the low state. These results give an interesting suggestion 
on the evolutionary state of the accretion disc in Circinus 
X-l. From the erratic and variable X-ray behaviour in the 
high state, the disc is considered to be in a later stage of 
the transition phase. The systems's X-ray behaviour is thus 
dominated by the wave induced accretion. In the low state, 
the disc is considered to be in the developing phase. Thus, 
the X-ray maxima resulting from the direct accretion are 
concentrated at periastron, as shown in Table |5] 



5 CONCLUSIONS 

In the framework of interaction between the Be star with the 
circumsteller disc and the neutron star, we have carried out 
the three dimensional SPH simulations of the long-term evo- 
lution of accretion disc around the neutron star in a copla- 
nar Be/X binary with a short period (P or b = 24.3, d) and 
a moderate eccentricity (e = 0.34), taking into account the 
phase-dependent mass transfer from the Be disc. Our main 
conclusions are summarized as follows: 

(1) A persistent accretion disc is formed around the neu- 
tron star in Be/X-ray binaries with a short period and a 
moderate eccentricity. 

(2) The disc has a one-armed structure of transient na- 
ture. The one-armed density wave is excited by the ram 
pressure of the material transferred from the Be disc at each 
periastron passage, and is damped by the next periastron 
passage. 

(3) The disc evolves via three distinct phases: the first 
developing phase, the second transition phase and the fi- 
nal quasi-steady state. In the developing phase, the disc is 
formed and develops towards a nearly Keplerian disc. Then, 
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in the transition phase, the approximately Keplerian disc 
grows with time. Finally, in the quasi-steady state, the mass- 
accretion rate on to the neutron star is on average balanced 
with the mass-transfer rate from the Be disc, and the disc 
exhibits regular orbital modulation. 

(4) There are two mechanisms that causes the orbital 
modulation of the accretion rate in Be/X-ray binaries: the 
direct accretion and the wave-induced accretion. While the 
former is due to the low specific angular momentum of the 
material transferred from the Be disc, the latter results from 
the one-armed wave induced by the periodically varying 
mass-transfer from the Be disc. The accretion rate has dou- 
ble peaks per orbit by these mechanisms in an earlier stage 
of disc evolution. In a later stage, however, it has only a 
peak arising from the wave-induced accretion. 

Be/X-ray binaries are an ideal laboratory of systems 
for studying physics of accretion. Unlike close binaries, most 
Be/X-ray binaries have eccentric orbits, the double circum- 
steller discs and can also be highly inclined systems. Such 
systems provide us a valuable opportunity to study the ef- 
fects of the periodically-changing tidal potential and mass- 
transfer rate and the inclination angle on the structure and 
evolution of the accretion flow, which hitherto have been 
studied little. In addition, the knowledge on these effects 
are applicable to a wide range of astrophysical objects, e.g., 
AGNs with periodic light curve such as OJ 287 which is con- 
sidered to have binary black holes with double circumblack- 
hole discs in the central region of nuclei. Much more work is 
desirable both theoretically and observationally in order to 
understand this interesting and important group of objects. 
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